clear all
set more off


local J = `1'

local N = `2'


use flexlogit



* prepare sample according to VU and z
clogit chosen x* z,group(ID)
gen vu=0
forvalues i=1/14{
replace vu=vu+_b[x`i']*x`i'
}
* find good with highest VU: good m
bys ID:egen maxvu=max(vu)
gen vuismax=(maxvu==vu)
gen rn=runiform() if vuismax==1
bysort ID:egen maxrn=max(rn)
gen vuismaxrn=(rn==maxrn)
drop vuismax
rename vuismaxrn vuismax
drop rn maxrn
gen vumaxID=altID if vuismax==1 
bys ID: egen vumaxmaxID = max(vumaxID)
drop vumaxID
rename vumaxmaxID vumaxID 

* find good with highest z: good 1
bys ID: egen maxz = max(z)
gen ismax = (maxz == z)
gen rn=runiform() if ismax==1
bysort ID:egen maxrn=max(rn)
gen ismaxrn=(rn==maxrn)
drop ismax
rename ismaxrn ismax
drop rn maxrn
gen maxID = altID if ismax == 1
bys ID: egen maxmaxID = max(maxID)
drop maxID
rename maxmaxID maxID 

* drop the choice set with the same good with highest VU and z
drop if maxID==vumaxID
* adjust terms, coefficients and ratio formula 




clogit chosen x* z, group(ID)

local naivelogit=_b[z]/_b[x1]



*Create attributes of rival goods
preserve
keep ID altID x1 x2 x3 x4 x5 z

forvalues i=1/5{
rename x`i' x`i'_
}

reshape wide x* z, i(ID) j(altID)
tempfile attributes
save `attributes', replace
restore

merge m:1 ID using `attributes'
assert _merge == 3
drop _merge




gen vijnaive=0
forvalues i=1/14{
 replace vijnaive=vijnaive+_b[x`i']*x`i'
}

replace vijnaive = vijnaive+_b[z]*z
*Create estimate of dv1dz1

gen dv1dz1est = _b[z]

local vij = "vijnaive"
local label = "naive"
do programs/makelogitprob.do "ID" "vij" "`label'" "`J'" "altID"

*Adjust weights so that denominator is never too small causing weights to blow up
replace sijnaive = max(sijnaive,.1)
replace sijnaive = min(sijnaive,0.9)
gen vusijnaive=0 
forvalues k = 1/`J' {
gen si`k'naive = sijnaive if altID == `k'
bys ID: egen si`k'maxnaive = max(si`k'naive)
drop si`k'naive
rename si`k'maxnaive si`k'naive
replace vusijnaive=si`k'naive if `k'==vumaxID
}
gen vubx=0
gen vubz=0
gen vub2x=0
gen vub2z=0
forvalues k = 1/`J' {



gen bx_`k' = si`k'naive/(1-sijnaive)

gen bz_`k' = si`k'naive/(1-sijnaive)

gen b2x_`k' = dv1dz1est*(1-2*sijnaive)*si`k'naive/((1-sijnaive)*(1+dv1dz1est*(1-2*sijnaive)*z))

gen b2z_`k' = dv1dz1est*(1-2*sijnaive)*si`k'naive/((1-sijnaive)*(1+dv1dz1est*(1-2*sijnaive)*z))

gen b2test_`k' = dv1dz1est*(1-2*sijnaive)*si`k'naive/(1-sijnaive)


*Adjust weights so that denominator is never too small causing weights to blow up

gen newval_`k' = 1+dv1dz1est*(1-2*sijnaive)*z
gen absnew_`k' = abs(newval_`k')
replace b2x_`k' = b2test_`k' if absnew_`k' <= 1
replace b2z_`k' = b2test_`k' if absnew_`k' <= 1
drop newval_`k' absnew_`k'


}

forvalues k=1/`J'{
	
replace vubx=bx_`k' if `k'==vumaxID
replace vubz=bz_`k' if `k'==vumaxID
replace vub2x=b2x_`k' if `k'==vumaxID
replace vub2z=b2z_`k' if `k'==vumaxID
}
sum bz_* bx_* b2x_* b2z_* 


do programs/makevariablessym.do `J'




clogit chosen x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 z zzmax zsum zzsum zxsum1 zxsum2 zxsum3 zxsum4 zxsum5 xsum1 xsum2 xsum3 xsum4 xsum5 vuz vuzz vuzx1 vuzx2 vuzx3 vuzx4 vuzx5  vux1 vux2 vux3 vux4 vux5, group(ID)


	
gen gammax = _b[vux1]
gen deltax = _b[vuzx1]


gen gammaz = _b[vuz]
gen deltaz = _b[vuzz]


local num = -_b[z]+gammaz+deltaz
local denom = -_b[x1]+gammax+deltax



****************************************FOD
*ratio of FOD at each choice set, ds1/dzJ // ds1/dxJ

gen numfod1 = -_b[z]*sijnaive*vusijnaive+(gammaz*vubz+deltaz*vub2z*z)*sijnaive*(1-sijnaive)
replace numfod1=. if ismax~=1
sum numfod1, meanonly
local numfod=r(mean)

gen denomfod = -_b[x1]*sijnaive*vusijnaive+(gammax*vubx+deltax*vub2x*z)*sijnaive*(1-sijnaive)
replace denomfod=. if ismax~=1
sum denomfod, meanonly
local denomfod=r(mean)

local theoryrat = `numfod'/`denomfod'



keep if _n == 1
keep chosen


gen naivelogit = `naivelogit'

gen theoryrat=`theoryrat'



drop chosen



